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ABSTRACT 


The  decay  of  earth  satellite b  has  been  analyzed 
using  both  a  numerical  integration  and  an  approximate 
method.  The  numerical  integretion  utilizes  the  vari- 
ation-of -parameter  method  with  second  order  earth 
gravitational  potential.  The  approximate  method  was 
developed  via  conservation  of  energy  equation  but  for 
unexplained  reasons  would  not  agree  with  observed 
satellite  decay  without  the  inclusion  of  an  additional 
linear  eccentricity  factor.  Including  (c/l00)®‘*‘‘  per¬ 
mitted  the  approximate  method  to  agree  with  numeri¬ 
cal  calculations  and  observed  data.  The  usual  re¬ 
striction  of  approximate  methods  at  low  eccentrici¬ 
ties  was  removed  by  the  addition  of  a  non-linear 
term  previously  computed  numerically.  The  linear 
approximate  method  was  found  to  be  valid  for  linear 
eccentricities  l.argcr  than  150,  The  effect  of  a  vari¬ 
able  drag  coefficient  and  projected  area  was  found  to 
bo  an  important  factor  and  of  the  order  e*  , 


TABLE  OF  CONTENTS 


Chapter 


ABSTRACT 


Ui 


LIST  OF  SYMBOLS 


vii 


I.  INTRODUCTION 


1 


II.  FORCES  ACTING  ON  AN  EARTH  SATELLITE  ....  2 


A.  Forces  due  to  the  presence  of 

celestial  bodies . 2 

B.  Forces  due  to  the  presence  of  the 

earth's  atmosphere.  . .  7 

III.  EQUATIONS  OF  MOTION  FOR  A  NEAR-EARTH 

SATELLITE . 2B 


IV..  NUMERICAL  SOLUTION 


33 


V.  AN  APPROXIMATE  METHOD  OF  SATELLITE  DECAY 

PREDICTIONS . 37 


VI.  DhSCUHSION  OF  SATELLITE  DECAY . '>3 

VII,  CONCLUSIONS . 6f> 


REFERENCES 


(.7 


LIST  OF  ILLUSTRATIONS 

Figure 

1  Relations  belsvccn  geocentric  and  orbit  axes  sysltMiiK  .  ,  4 

Z  Atmospheric  density . V 

3  Most  probable  molecular  speed  aa  a  function 

of  altitude . 13 

4  Minimum  molecular  speed  ratio  as  a  function 

of  altitude . 14 

5  Relation  of  re-emission  speed  ratio  to  incidtnt  speed 

ratio  yn  an  insulated  cylinder . 16 

6a  Trajectory  coordinates . 31 

6b  ];istantnneous  ellipse  coordinates . 31 

7  Specular  drag  coefficient  for  a  sphere  as  a  function 

of  al  ti  lude . 34 


TR  46S 


V 


LIST  OF.  ILLUSTRATIONS  -  (Continued) 

Figure  Page 

8  Density  scale  height  for  various  altitudes  ••••••  38 

9  Ratio  of  change  in  linear  eccentricity  to  change  in 

perigee  as  a  function  of  linear  eccentricity . 40 

10  Satellite  decay  density  parameter  . . 42 

11  Lifetime  of  an  earth  satellite  in  a  circular  orbit  ....  43 

12  Straiglu  line  approximation  to  linear  eccentricity 

vs  N.j.(Cjy\/w)  . . .  44 

13  Satellite  lifetime  prediction  as  a  function  of  altitude 

and  linear  eccentricity  . . 47 

14  Comparison  of  observed  and  predicted  s.atellite  lifetime 

for  various  altitudes . 48 

15  Norinali/iCd  altitude  factor . 49 

16  Comparison  of  observed  and  predicted  satellite  lifctitnc 

for  various  Cj^A/w . 50 

17  Decrease  of  linear  eccentricity,  . . 54 

18  Motion  of  argument  of  perigee  as  a  function  of  jiolar 

angle . .  .  57 

19  Motion  of  true  anomoly  at  low  altitudes  . . 58 

20  Motion  of  true  anamoly  at  low  altitudes  . . 59 

21  Variation  of  eccentricity  atdem'so . 60 

22  Decrease  ii\  perigee  at  demise . 61 

23  Variation  in  altitude  at  demise . 62 

24  Satellite  polar  angle  versus  altitude  ........  63 

25  Error  in  linear  eccentricity  encountered  by 

Cq  =  2.0  as  compared  to  cylinder  . . .  64 

LIST  OF  TABLES 

Tables 

1  Temperature-Height  Profile  of  1959  ARDC 

Model  Atmosphere . .  . . 10 

5- 1  Ballistic  Coefficients  of  Several  Satellites . 52 

6- 1  Effect  of  Drag  on  the  Motion  of  Orbital  Elements  ....  56 

6-2  Minimum  values  of  Eccentricity . 65 

TR  465  Vi 


LIST  OF  SYMBOLS 


Satellite  angle  of  attack 
Semi-major  axis 
Satellite  angle  of  yaw 
Magnetic  field  strength 
Satellite  angle  of  roll 
Drag  coefficient 

Ballistic  coefficient,  cm* /Kg 

Ballistic  coefficient,  ft* /lb 
Linear  eccentricity 

Jeffrey's  coefficient  for  the  third  order  gravitational  term 
Energy  due  to  drag  force 


E  or  o 


I' 

f 

G 


K 

H 

n, 


lo 

h 

J 


j 

KM 

K 

n 

L 

m 

M 


m 

N 

N 

s 


Eccentricity 
Perturbing  force 
'I'rue  anomaly 
Density  scale  i»eight 
Local  gravity 
i’ressure  scale  height 
Momentum  vector 
Altitude 

Modified  Bessel  function  of  zero  order 
Modified  Bessel  function  of  first  order 

Jeffrey's  coefficient  for  the  second  order  gravitational  term 
Current 

l':arth'8  gravitational  constant 

Knudsen  number 

dT,  /dH 
M 

Air  mass 
Satellite  mass 
Revolution 

Total  number  of  revolutions  in  spiral  decay 


VI  i 


TR  465 


r 

S 

'I' 


M 


T 


Hkiii 


U 

V 
Vi; 

I' 

V 

q 

V 

V 

r 

a 

r. 

q 

0 

\ 

p 

‘t> 

u 

L> 


.a 

T 


Orbit  parameter 

Legendre  polynomial 

Anomalistic  period 

Apogee  distance  from  earth's  center 

Perigee  distance  from  earth's  center 

Earth's  equatorial  radius 

Accommodation  coefficient 
Universal  gas  constant 
Radius 

Molecular  speed  ratio 
Air  temperature 
Molecular  temperature 

Satellite  surface  temperature 

Gravitation  potential 
Satellite  velocity 

Satellite  velocity  in  circular  orbit 

Satellite  velocity  at  |Jcrigoe 

Most  probable  molecular  velocity 
Molecular  re-emission  velocity 

Orbit  inclination 

Prolate  ellipsoid  coordinate 

Prolate  ellipsoid  coordinate 

Angle  between  current  and  magnetic  lines 

Longitude 

Atmospheric  density 
Polar  angle 
Latitude 

Angular  velocity 
Argument  of  perigee 
Arguincnt  of  right  ascension 

Vernal  equinox  (Ram's  horns  or  the  first  point  of  Aries) 


T  K  4  6  a 


tJ  H 


Subscripts 


0  Reference  altitude 

Base  of  a  constant  L 

1 

Specular  reflection 

Diffuse  reflection 
Total 


<!>' 

a 

E 

N  ► 

.Oo 


Represent  vector  direction 


I.  INTRODUCTION 


The  motion  of  satellites  is  no  longer  a  new  subject.  The  flood  of 
papers  and  notes  in  current  journals  would  indicate  that  considerable  at¬ 
tention  has  been  given  to  the  problems  of  space  flight  and  atmospheric 
re-entry.  In  many  of  these  studies  the  analysis  pertains  to  one  of  two 
aspects:  orbital  motion  or  re-entry.  Less  attention  has  been  given  to 
discussion  of  the  satellite's  motion  throughout  its  entire  lifetime.  In 
jjarticular,  analysis  of  the  transition  flight  bt-twetm  orbit  and  re-entry 
has  been  lacking.  '1‘his  study  has  been  a  two-prongecl  analysis  of  satel¬ 
lite  iiiotioii  witli  tile  object  being  decay  prediction.  Tiie  analy.sis  is  two¬ 
pronged  in  that  lioth  an  iipproxini.itn  metliod  which  permits  gr.ipliieal  or 
simil.irly  c([ually  e.isy  solutions  and  refined  methods  recpiiring  inimeriiail 
integration  Ijy  modern  (;omputing  macliines  are  employed, 

'I'lu!  .approximate  nietliud  reviewed  .anti  improved  liy  an  cmpiric.al 
f.aclor,  is  liiniteti  in  long  r.ange  pretiiction  liy  tii(‘  n.aturt!  of  tiie  sun's  in¬ 
fluence  upon  tiie  e.artli's  atmosphere.  Thi.s  is  importanl,  liowtwer,  only 
.aiiove  iOO  km.  The  form  t>f  the  solutions  obtaineti  liiianigii  tiie  apjiroxi- 
mate  nietiiod  is  well  suiteti  for  analysis  of  tiie  gro.ss  snlelliti?  inotion,  tie- 
cay,  anti  atmuspiieric  density.  Mt»re  tiet.ailed  infornialion  is  to  iie 
giatintretl  from  tiie  tiunierituil  solutions  of  tiie  gener.ai  tli ree -dimeii sioiiai 
(•((uations.  'J'iie  ttiree-tiiiiiensitjiial  etpintions  of  motion  etpiate  tiie  s.atel- 
lite  acceleration  to  the  perturljiiig  function,  Tiie  pe rturiiing  fumitiun  is 
m.acie  up  of  tiie  first  and  second  orders  of  tiie  earlli' s  gravitational  poten¬ 
tial,  and  lift  and  drag  forces.  Tlic  orientation  of  tiie  satellite  in  sjiace  is 
also  noteti.  Tiie  set.-ond  order  differenti.a]  etpiatitins  are  reduced  to  six 
f'ir.st  tirder  ilifferential  etiuntioiis  using  file  well  known  variation -of - 
parameter  technique.  The  solution  of  these  six  equations  siinullaneously 
is  ol3taiiie<i  numerically  with  an  IBM  709.  Tiie  viiri.ation-of -paranictc r 
perturbation  tcclinique  is  w<-!l  suited  for  describing  the  satellite  motion 
througliout  tin-  tr.insition  fliglit  but  is  limited  lo  the  order  of  several  hundred 
rovoliitinns  by  madiine  time.  Sinc<3  this  is  a  numerical  solution,  the 
.inalysis  of  satellite  motion  consists  in  observations  of  several  cases,  This 
c.ilculation  desi  rilies  tile  secol.ar  and  pt-Tiodic  motion  of  the  orbital  elements. 
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The  numerical  solution  also  serves  as  an  "exact"  solution  to  check  the 
approximate  method.  Also  the  approximate  method  is  corrected  to 
agree  with  the  numerical  solution. 


II.  FORCES  ACTING  ON  AN  EARTH  SATELLITE 


A  satellite  launched  into  an  orbit  about  the  earth  travels  through 
a  magnetic,  field,  a  radiation  fi<'ld,  and  the  earth's  atino.sphere.  These 
fields  can  cause  either  a  retarding  or  an  accelerating  force.  The  most 
important  force  is  the  drag  and  lift  created  as  the  satellite  passes  through 
the  earth's  den.se  atmosj)hcre.  Each  of  these  forces  will  be  examined 
in  d(!tail  in  this  section.  Four  fundamental  satellite  shapes  will  be  con¬ 
sidered;  spher<!,  flat  |)late,  cyliiuler,  and  prolate  ellipsoid.  Most  any 
actual  satellite  ciiu  be  tluuight  as  a  composite  of  one  or  more  of  these 
shaiHiH.  The  final  outcome  of  this  section  is  the  writing  of  the  iierturbing 
function, 


A.  Forces  (Uu!  to  tlu;  presence  of  celestial  bodies 
1,  'I'he  (‘artb 

Tbt!  gravitational  potential  of  the  tnirth  can  i)0  expressed 
as  an  infinite  series 


r 


1  - 


00 


r 


n=:l 


n 

P  (sinv'/) 
n 


(i.l) 


where  (//  represents  the  latitude  and  R  the  earth's  equatorial  radius.  The 
equatorial  radius  is  given  as  6378.145  km.  The  earth's  gravitational  con¬ 
stant,  KM,  is  equal  to  3.98614  X  10®  km®  /  sec^  .  The  coefficients  of  the 
series  were  originally  designated  by  Jeffrey  (1)  as  .1  and  D,  where 
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The  J  and  D  coefficients  have  been  redetermined  since  that  time  by 
Cornford,  King-Helc,  and  Merson  (Ref,^)  through  observations  of  the 
secular  satellite}  orbit  motions.  The  rurrenl  values  are 


J  =  ( lf),t4.6  i  0.3]  X  lO'*^ 

D  =  (5.7  1  0.8]  X  10‘‘ 

Itefore  the  secular  motions  of  the  earth  satoJlito  were  available 
for  stufly,  the  odd  h.armonic  term.s  of  tlie  series  for  the  earth's  gravita¬ 
tional  potential  were  assumed  to  be  zero,  i.e,,  the  earth  was  thought  to 
be  Hynim«!tric  ai)Oul  t)ie  equator.  Now  it  has  been  sliown  tli;it  the  odd 
harmonics  do  exist,  'fhe  first  harmonit;  can,  however,  be  tna<!o  zero 
by  plat  Ing  the  axis  system  at  thtr  t:cnter  of  mass.  When  tht;  center  of 
mass  is  not  at  the  point  of  symmetry,  other  diffiiuiltit's  will  be  en- 
counttu-ed  in  determining  the  height  frojii  the  eartli's  surface  and  the  cor- 
rt^sponding  dtuisity. 

Tlie  force  per  tuiit  mass  acting  on  the  satellite  due  to  the  earth's 
gravitational  potential  is  the  potential  gradient  T  -  •  VU.  Willi  spherical 
polar  cooiulinate s  cenlereti  .at  the  center  of  m.iss,  thi;  components  of  the 
gravitational  force  along  the  radius,  latitude  an<i  longitude,  respectively, 
a  re 


F 

r 

"nT 


ou 

Or 


(1 


)  ... 


—  30  sin*  Ip  +  3)  +  . 


(2.2) 


^iP 

m  r  dtf 


=  JKM 


sinZ^  + 


sin  (3  —  7  sin*  ip  )  + 
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m  “  r  sin^ 

It  Is  usually  more  convenient  to  use  a  coordinate  system  located 
in  the  orbit  plane  rather  than  the  geocentric  system  above.  The  axis 
system  in  the  orbital  plane  is  designated  r,  ()>,  a  (see  Figure  1).  The 
direction  of  (|)  is  in  the  plane  perpendicular  to  the  radius  vector  and  a 
is  the  direction  normal  to  the  orbital  plane.  Several  relations  between 
the  two  axis  systems  will  be  useful.  These  can  be  found  through  the 
use  of  spherical  and  plane  trigonometry. 


N 


Figure  1,  Relations  between  geocentric  and  orbit  axis  system 
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These  relations  are 


( I )  cos  4>  —  cos  (X.  —  n)  cos  ^ 
(1)  sin^  =  sin4>  sina 
(3)  COSO  =  cos(^  sin 


(4)  8in(J 


8in(X.  —  n) 


I  tani> 

(5)  tana  = 

With  the  axis  system  located  in  the  orbital  plane,  the  components  of  force 
du(!  to  the  Ccarfii's  gravitational  potential  are 


m  r‘ 


30sin*  (f  a  +  3 )  + 


Tn'  ~  ^  ”  ■*‘**”*'1’  (^)  sin*  a 

—  30sin*  (f  a  +  3)  +  •  •  •  j 

— .i  =  ,dj  KM  sin^a  sin  (j^cos <)»  +  ^  KM  j^<18sin*a  sin*  <|)  cos  (j)  (2.3) 

—  12 sin*  a  sin  (|i  cos  (j>  j 

^  a  H*  D  R*  r 

TfT  ~  J  -pc  sin<t)  siao"  cos  a  +  -y-  -pt-  KM  |  28  sin*  <|)  .sin*  a  cos  a 


128in<|)  sin  a  cos  a 


TR  465 


5 


In  addition  to  ilic  forces  due  to  the  earth's  gravitational  potential, 
there  is  a  force  due  to  the  interaction  between  the  satellite  and  the  earth's 
magnetic  field.  This  force  was,  for  example,  predominant  in  causing  the 
motion  of  the  Tiros  I  spin  vector  (Hef.3).  Previously  (Ref.  4),  the  only 
interaction  considered  was  due  to  the  cutting  of  the  earth's  magnetic 
lines  by  the  inducetl  electrical  charge.  Another  interaction  which  should 
bo  included  is  between  the  satellite  internal  current  and  the  earth's  mag¬ 
netic  field.  Such  a  force  is  described  by 

F  =  j  X  n, 

In  particular 


dF  =  0.10i(dSX  13,  )  dynes 

where  the  current  (anipercs),  i,  flows  in  .a  leiigth  (centimeters)  dS  of 
wire  in  a  magnetic  field  (Gauss)!!.  For  Tiros  1,  there  is  an  equivalent 
curr<!nt  of  om?  ainpcire  flowing  around  the  eight-foot  base  diameter  which 
proiiuctfs  a  force  of 

F  =  Zi  sine  0  dynes 

where  0  is  th(.'  .nnglc  betwt  en  tin;  current  ami  the  magnetic  lines.  This 
fnin:e  is  compa  rabi <!  to  the  dr.ig  force.  The  force  is,  however,  not  always 
of  the  same  sign.  The  existi'nce  of  this  force?  resulting  from  an  internal 
current  is  extremely  important  and  bears  further  examination. 

Z,  The  Sun  and  Moon 

The  sun  exerts  considerable  influerce  upon  a  satellite, 
cither  by  gravity,  radiation  pressure,  or  indirectly  through  its  influence 
upon  the  earth's  atmosphere.  The  nature  of  these  influences  has  been 
discussed  in  detail  by  many  authors.  The  conclusions  to  date  are  that 
atmospheric  heating  is  the  means  of  increasing  the  density  during  exposure 
to  the  svm.  Another  observed  phenomenon  associates  an  increased  drag 
with  solar  storms.  This  latter  correlation  has  lessened  the  hope  of  de¬ 
veloping  long  range  decay  prediction  since  solar  storms  are  as  yet  un- 
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predictable.  The  obvious  escape  from  this  difficulty  is  to  limit  the  analy¬ 
sis  to  altitudes  below  300  km  where  the  sun's  i:ifluence  has  been  ob¬ 
served  to  be  small  or  non-existent.  Other  difficulties  are  encountered 
-also  above  300  km.  At  extreme  ranges,  of  the  order  of  the  earth's  radii 
or  more,  the  solar  gravity  field  and  radiation  pressure,  as  well  as  the 
lunar  gravitational  field,  must  be  included  in  the  perturbing  function. 

The  sun  and  moon  have  an  influence  on  the  earth's  at- 
mosphere  below  ^00  km  in  the  form  of  atmospheric  tides;  Atmospheric 
tides  have  been  detected  from  pressuremeasurements  to  occur  twice  a 
day  ami  have  ajr~amplitude  of  2.  mm  llg  prerssure  change.  This  is  equiva¬ 
lent  to  a  li  km  altituth'  changOj^  Obsi^yationj^of  satellites,  Jiavc  not  as  yet 
jr^evealed  a  liehavior  which  c.'in  l)e  attributed  to  atmospheric  tides.  --Tt  is-un 
likely  tli.it  a  variation  of  2.  km  would  be  important  in  view  of  the  present 
stall!  of  knowledge  of  the  varying  ionospheric  lieusity. 

il.  Forces  due  to  the  presence  of  tlic  earth's  atmosphere 

The  satellite  passing  through  the  eartli's  .itmosphurc  experiences 
aerodynamic  forces  of  lift  and  drag  as  v/ell  as  moments  about  the  center 
of  gravity.  These  acruiiyiuimic  forces  and  moments  can  seriously  effect 
llui  satellite  orliil.  'I’he  effect  of  drag  is  well  known  but  whether  lift  and 
moments  augment  or  diminish  the  decay  is  not  so  well  understood.  An 
analysis  of  the  effects  of  these  three  qu.intities  on  the  lifetime  of  a  satel¬ 
lite  is  presented  in  this  section.  First,  however,  the  earth's  atmosphere 
should  be  specified. 

In  free  molecular  flow,  the  number  and  energy  slate  of  the  air 
particles  is  of  interest.  Usually,  the  density  is  expressed  functionally 
as  an  exponential 

p  =  Po  exp  [-  (2.5) 

where  H  represents  the  scale  height  defined  as 
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„  RpT 

”  “  gM  -  "iWT 

As  can  be  seen,  the  scale  height  is  a  function  of  the  atmosphci'ic  tempera¬ 
ture  and  molecular  mass.  “  Consequently,  Equation  (d. 5)  is  limitcd_to  an 
altitude  close^lo  the  reference  altitude  (hp  ).  Above  300  km,  the  variation 
of  the  scale  height  and  density  are  nearly  linear  and  Equation-(.l ,5)  is  - 
valid  over  a  range  of  100  km  above  the  reference  altitude.  As  an  .example,, 
starting  at  300  km,  tho  density  using  Equation  (,i.6)  will  be  ,10  per  cent  too 
low  at  360  km.  The  error  is  much  larger  if  the  reference  altitude  is  be¬ 
tween  100  and  300  km.  'I'iiis  is  to  be  expectoei  since  the  region  between 
100  and  300  km  is  the  transition  between  two  functions  which  fit  the  at- 
mospheric  data.  As.a  fur.tlier  .example, _for  a  reference  altitude.at  .,100  km, 
the  error  at  ,160  ktn  is  about  .‘>0  perTent  too  low. 

An  alteiMiate  expression  wliich  more  closely  fits  observed  values 
above  3  50  km  is  .a  powitr  function. 


The  value  of  n  whi<:h  fits  the  1959  AKDC  model  atmosphere  is 
about  1/6.  A  smoother  r«!pro8ont;,tion  would  result  if  the  value  of  6.38 
i.*,  uae*!.  Figure  f.  is  a  comp,arison  of  this  fuiuTiun  with  the  atmos|)heric 
mudel, 

In  place  of  eitlier  of  these  two  equations,  tlie  actu.aJ  atmospheric 
model  can  be  used.  In  machine  computations  such  as  used  for  Ibis  study, 
•I  subroutine  to  look  up  and  interpolate  the  atmospheric  table  can  be  in¬ 
corporated.  This  form  is  necessary  in  attempting  to  preserve  accuracy 
throughout  the  compulations. 

Each  of  these  three  methods  inlicrently  assumes  that  the  earth's 
atmosphere  has  the  same  shape  as  the  earth.  Logically,  it  would  be 
expected  that  the  oblateness  of  the  atmosphere  would  be  different  from  the 
oblatcnesB  of  the  earth.  A  more  exact  shape  is  currently  being  evaluated 
by  Jacchia  (Ref,  5). 

In  addition  to  the  density,  a  model  of  the  atmospheric  tempera¬ 
ture  and  mass  is  required.  These  two  quantities  are  sometimes  grouped 
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together  as  the  molecular-scale  temperature.  The  following  is  the 
general  form  of  the  molecular-scale  temperature  function  (Ref.  6), 


M 


b 


where 


n  =  Cicopotcntial  altitude  in  km 

a  The  molecular-scale  temperature  is  °K  at 
altitude  II 

II|^  =  Goijpotential  in  km  at  the  base  of  a  particular 
layer  charactoriined  by  a  specific  value  of 


(  '  kv)  =  '1  he  value  of  T.,  at  altitude  II, 
M  M  b 


The  values  of  the  coefficients  arc  given  in  Table  1. 


(2.8) 


Table  1 

I'emperature -Meiglit  l^rofile  of  1959  ARDC  Model  Atmosphere 


!'•> 

<'^M>b 

1,..  -  “K/km 
M 

79 

1 65.66 

0.0 

90 

165.66 

4.0 

105 

ZZ5.66 

Z.O 

160 

13Z5.66 

1.0 

170 

14Z5.66 

0.5 

ZOO 

157  5.66 

0.35 

700 

33Z5.66 

0.35 

1 .  Drag  forces 


The  drag  force  acting  on  the  satellite  is  proportional  to  the 
dynamic  pressure  and  the  area.  The  factor  of  proportionality,  termed  the 
drag  coefficient,  is  a  function  of  the  relative  velocity  of  the  satellite  with 
respect  to  the  atmosphere,  the  density  of  the  atmosphere,  and  the  surface 
conditions  of  the  satellite.  The  density  of  the  atmosphere  dictates  which 
aerodynamic  approximation  is  used  to  describe  the  drag  coefficient.  The 
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regions  of  flow  to  which  each  approximation  is  applicable  are  commonly 
separated  by  the  Knudsen  number,  a  non-dimensional  ratio  of  the  mean 
free  path  to  some  pertinent  body  dimension,  Experiments  in  free 
molecular  flow  (Ref.  7)  indicate  that  the  "so-called"  free  molecular  theory 
becomes  applicable  when  the  Knud-en  number  is  greater  than  two.  As 
discussed  previously  (Ref,  4),  free  molecular  theory  is,  on  this  basis, 
valid  for  all  satellite  orbits  above  130  km  (assuming  a  satellite  length  of 
is  m).  The  launching  of  larger  satellites  such  as  the  100-foot  balloon, 
will  require  the  addition  of  a  transition  drag  correction  to  the  free 
molecular  flow  theory  as  the  s.atcllites  descend  below  160  kni.  The  de- 
pcn<ionce  of  the  drag  coefficient  on  the  satellite  velocity  and  surface  con¬ 
dition  will  not  be  derived  from  first  principles  in  the  succeeding  sections. 
For  a  more  fletailed  discitssioti  of  the  dr.ag  coefficient,  the  reader  is 
ref (jr rod  to  Refs.  8,  9,  and  10, 

(a)  Estim.'itiun  of  molecular  spi-ed  ratio 

One  of  th«‘  ba.sic  dimensionles.H  parameters  in  the  free 
inolecuhar  flow  is  (■■•illc*!  th<!  molecular  spcjed  ratio,  and  is  defined  ,i8! 


v 


whtM'e 

'I'lie  satellite  velocity  is  denoted  as  V  and  the  most  probable 
molecular  speed  as  v.  The  inolecula r- scale  temperature,  and  the  gas 
constant  R  are  obtained  from  Rt;f,  6. 

A  range  for  possible  molecular  speed  ratios  encountered 
by  satellites  at  various  altitvides  is  of  interest.  The  minimum  molecular 
speed  ratio  will  be  defined  in  terms  of 

Vc 

S  .  =  -  — 

min  V 

where  Vc^  is  the  circular  satellite  velocity  at  a  distance  r  from  the  earth's 
center  and  is  defined  as 


1  1 
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given  as 


For  vehicles  in  elliptical  orbits,  the  total  velocity  is 


■  y  _  /km  (1  +_2e  cosf  +  e^ )  ^ 
“  y  -r  'T'+  e  cbsl 


_Tho_earth.'s  gravitational  constant  is  represented  as  KM.  The  velocity 
at  perigee  can  be  expressed  as  a  function  of  the  circular  velocity  by 


V  =  Vc  (1  +  e)  (2.10) 

q  - r-  -  - . -  —  ■ 

It  can  be  seen  that  the  velocity  at  perigee  for  an  elliptical 
orbit  of  low  eccentricity  (such  as  most  satellites  around  the  earth)  is 
only  slightly  larger  than  the  circular  satellite  velocity  at  the  sajne  point, 
and  that  the  minimum  8p<!ud  ratio  is  a  good  representation  of  the  actual 
molecular  speed  ratio  at  the  perigee  points. 

The  molecular  velocity  as  a  function  of  altitude  is  il¬ 
lustrated  by  Figure  3  while  the  values  of  as  functions  of  altitude 

are  plotted  in  Figure  4.  It  is  seen  that  in  the  150-600  kin  altitude  range, 
the  range  of  interest  for  present  earth  satellites,  the  minimum  molecu¬ 
lar  speed  ratio  ranges  between  six  and  ten.  The  most  probable  molecu¬ 
lar  speed  increases  from  a  value  of  0.782  km/soc  at  150  km  to  1.285  km/s 
at  600  km. 

For  conditions  of  diffuse  reflection,  the  drag  coefficient 
will  be  based  on  the  re-emisiion  speed  ratio,  S^,  The  velocity  of  re- 
emission  is  not  clearly  defined.  Most  likely  it  is  a  function  of  the  sur¬ 
face  conditions.  However,  the  only  readily  available  measure  of  condi¬ 
tions  is  the  surface  temperature.  Consequently,  it  is  usual  to  assume  the 
velocity  of  re-emission  to  be 
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ALTITUDE  (K 


®  ®  7  8  9  10 

MINIMUM  MOLECULAR  SPEED  RATIO,  Sm,n 
Figure  4.  Minimum  molecular  speed  ratio  as  a  function  of  altitude 
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The  wall  temperature  :varieb  from  recovery  temperature 
to  ambient  gas  temperature.  Above  an  altitude  of  about  200  km  the  heat 
tranafer  due  to  the  incident  molecules  becomes  less  than  the  solar  flux. 
Consequently,  the  wall  temperature  is  detozmiinod  from  a  radiation 
balance.  Since  it  is  unlikely  that  the  emissivity  of  the  satellite  can  be 
made  less  than  0,01,  the  maxintum^.surface  tc*mperaturu  above-2 00  km-  ^ 
will  bo  1250  "K.  This  temperature  will  ahvay.s  be'lesH  than  the  free- 
stream-temperature  above  200  km  in  which  case  the  re>emi8sion  speed 
ratio  can  be  speenfied  as  equal  to  the  incident  mplecular  speed  ratio. 

!•  or  flight  below  200  km  whore  free  molecular  heat  trans¬ 
fer  is  important,  the  re-omission  speed  ratio  c.nn  be  computed  (Kef,  7) 
for. the- extreme  case  of  no  heat  transfer.  This  neglects  effeets-of  radi¬ 
ation  and  eon<liiition.  Staldcr  obtained  experimental  veriflcatioii  for  this 
condition  UHing  a  Knudsen  number  of  two.  The  functional  relalioii  hi'tw(,UMi 
the  two  Mpeofl  vatios  can  bo  approximjited  as 


Sj.  «  1 .64  5  •  .  0.18(6  -  .S)2 


S  <  6 


(^..11) 


This  is  to  say,  the  re-emission  speed  ratio  is  always  much  smaller  than 
the  incident  molecular  speed  ratio  (see  Figure  5).  The  effects  of  conduc¬ 
tion  and  radiation  are  expected  (o  dra.stically  reduce  the  surfai c  tempera¬ 
ture  and  thereby  increase  the  re-emis.s.ion  speed  ratio.  If  recovery  tem¬ 
perature  was  reached  or  in  effect,  radiation  and  conduction  were  negligible, 
the  drag  due  to  the  re-emission  could  conlrilmte  a  maximum  of  4U  per  cent 
of  the  total  drag  force  as  will  be  seen  shortly.  In  practice  this  will  not 
occur  but  it  is  interesting  to  note  that  the  drag  coefficient  could  become 
atmospheric  dependent  through  a  variation  of  the  re-emission  speed  ratio. 

(b)  Estimation  of  drag  coefficient 

Initially  the  satellite  will  be  in  free  molecular  flow  which 
is  assumed  to  exist  when  the  atmospheric  mean  free  path  is  equal  to  twice 
the  characteristic  dimension  of  the  satellite.  Experiments  have  been  made 
which  verify  this  assumption.  When  the  Knudsen  number  becomes  loss 
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than  two,  some  type  of  transition  flow  such  as  slip  flow  occurs  until  con¬ 
tinuous  flow  is  reached  at  Knudscn  numbers  of  about  K  “  0.0<J  for  or- 

n 

bital  vehicles.  This  transition  region  will  most  likely  begin  as  the  satel¬ 
lite  is  still  orbiting.  As  an  example,  the  transition  point  (K^  =  <J.O)  for 
the  Discoverer  satellites  is  at  abo'*t  130  km.  Below  130  km  drag  coef¬ 
ficients  based  on  experiment  or  one  of  the  theories  such  as  those  which 
include  one  collision  between  air  molecules  would  be  used.  In  any  case,  the 
general  trend  is  to  decrease  the  drag  coefficient  from  the  free  molecular 
condition.  As  will  be  mentioned  later,  at  an  altitude  of  130  km,  the  satel¬ 
lite  is  about  to  begin  its  terminal  phase.  Numerical  calculation  using  free 
molecular  flow  drag  coefficients  indicates  the  satellite  traverses  3/4  or 
loss  of  an  orbit  in  descending  from  130  to  90  km  (see  P'igure  <i4).  The  er¬ 
ror  involved  in  using  free  molecular  drag  coefficients  in  place  of  the  more 
realistic  slijj  or  continuum  flow  could  not  cause  an  error  greater  than  two 
revolutions.  This  is  usually  negligible  for  lifetime  predictions.  On  the 
other  Ijantl,  if  an  impact  point  is  desired,  the  more  exact  drag  coefficients 
will  have  to  bo  used. 

The  basic  assumption  of  free  molecular  flow  is  that  t)iu  air  moli.'cules 
liave  a  Maxwellian  velocity  distribution  superimposed  upon  the  uniform  mass 
vehjcity,  Then,  by  n«'glecting  collisions  l»ctween  molecules,  the  flux  in- 
citlence  upon  a  given  area  is  compul<‘ci,  A  firag  forc;(!  also  occurs  in  the  re¬ 
flection  of  molecules  at  tJu;  surface.  It  is  n(!ccssary,  however,  to  make  an 
estimate  of  the  surface  conditions  in  order  to  i>reciict  tlie  mode  of  reflection. 
Usually  a  combination  of  specular  ;ind  diffuse  reflection  occurs.  'I'he  drag 
ctxTficients  tor  four  commonly  eiicountcTcd  .shapes  are  {lescribed  below. 

i.  .Sphere 

With  specular  reflection  from  the  surftace,  the  free  molecular 
drag  coefficient  for  a  spiicrcr  based  upon  the  frontal  area  is 


C 


D 
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Fox'  a  diffuse  reflecting  sphere 


°ti  ^/l;  s‘ 


2*/ff 

"JS“ 


(2,13) 


here  is  the  re-omittx.-d  speed  ratio, 
'.('he  coiTibincd  drag  coefficient  is 


8  r 


(-i.H) 


wher(!  is  tlic  proportionality  const.aiit  hc-tweiui  diffiioe  and  specular 
r idlections .  In  general,  this  coefficient  is  assumed  to  have  a  value  hc- 
twexui  0.8  and  1 .0. 

ii.  Flat  plate 

'I'hongh  no  siitellites  presently  in  orbit  arc  approximated 
as  flat  plates,  neverthele.ss  certain  portions  of  tht?  siatullitcs  arc  fl.at. 

As  an  cxainph;,  the  tumbling  cyliinler  will  at  frirtiueiit  intcrv.als  project 
a  flat  (uul  to  the  air  stream,^ 

'I'lie  drag  (  oefficixnits  for  a  fl.it  jxlale  at  an  angle  of  attack, 
a,  with  tlu;  wind  velocity  ''cetor  is 


Cl,  ~  sin  aJ  —  sin  A  exjxi  --  S*  sin^  A  )  +  4  (  sin*  A  + 
‘  s  L  TT  S  \ 


+  2^^  (S  sin  A)  j 


(2.15) 


Cj^  =  — ~ —  exp  (  —  S*sin  A  )+2sinAfl+  )  erf(S  sin  A) 
d  Kfli  S  '  ' 


(2.16) 


+  sin*  A 
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The  total  drag  coefficirnt  is 


a  s 

iii.  Cylijidor  transverse  to  the  flow 

The  moat  importiint  basic  shape  which  closely  resembles 
niany  presently  orbiting  satellites  is  the  cylinder. 

The  drag  coefficient  resulting  from  flow  transverse  to  the 
cylinder  based  on  tlie  frontal  area  is  described  below. 


.  (1  + 

+  — - - 


I-  I, 


55"" 

r 


isT ff 

<?xp 


I-  (1  i  .18^1, 


(a.lH) 


(ri.19) 


A.S  in  all  cases,  the  toUil  drag  is 

‘  il  s 

iv,  l^rolale  ellipsoid 

'I'his  .shape  resembJes  some  of  the  re-entry  conc'S  and  may 
be  useful  at  times.  The  surfaiM?  area  of  a  prolate  ellipsoid  is  cxprcsseti 
as 

f*  \  I  7 

A  =  J  C2  [(v^  -ti2)(^2  -  1)]  dr]d<\> 


TH  46  5 


19 


sin  0  cos  0  I  IV 


D 


specular 


The)  elruK  cocfficiontH  for  the  four  basic  shapes  have  lieun  presented 
in  terms  of  tlie  ir.olecular  speed  ratio  and  the  satellites  projuctetl  area. 

What  is  H^’neirally  done  is  to  linoari/.e  the  functional  representations  by 
asHumin)'  S  >  >  1.  'I'liis  corresponds  to  hypersonic  flow  and  nj>proache8 
the  Newtonian  theory.  This  assurnption  is  only  valid  if  liie  iiccuracy  do- 
Hiro<l  is  to  the  order  of  the  orbital  eccentricity.  For  hi(»her  order  ac¬ 
curacy,  such  as  that  of  or  higher,  it  is  essential  that  more  exact  repre¬ 
sentations  of  the  free  molecular  drag  coefficients  he  employed. 

The  form  of  tlic  drag  coefficient  is  only  a  part  of  tlje  difliculty 
encountered  in  analyzing  tlic  satellite  drag.  The  orientation  of  the  satel¬ 
lite  itself  is  also  important.  However,  when  the  hypersonic  flow  assump¬ 
tion  is  made,  the  drag  coefficient  per  unit  frontal  area  is  identical  for  the 
sphere,  cone,  and  cylinder.  Thus,  when  the  analysis  is  only  to  order  of  the 
eccentricity,  the  satellite  orientation  is  not  important.  In  the  attempt  at  re¬ 
fining  the  satellite  decay  predictions,  higher  order  terms  are  included  in  the 
orbit  equations  and  in  keeping  with  this  increased  precision,  the  satellite 
orientation  as  well  as  the  complete  expression  for  the  drag  coefficient  must 
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be  indicated.  In  the  following  paragraphs  attention  will  be  given  to 
evaluating  the  drag  coefficient  on  an  actual  satellite. 

The  most  evident  way  to  begin  is  to  .analyze  the  space  and  earth 
stabilized  satellites.  At  one  time  it  was  thought  that  a  satellite  could 
be  stabilized  in  sp.ace  in  much  the  same  manner  as  a  gyroscope.  This 
assumption  Nvas  not  upheld,  as  was  proved  by  the  observation  of  the 
rapid  motion  of  the  spin  axis  of  Explorer  IV  about  the  celestial  sphere. 
This  posed  a  dilemma  for  photographing^  the  earth  from  Tiros  I.  It  has 
been  siibstantially  demonstrated  (Ref.  3)  that  effects  of  differential 
gravity  forces  cannot  be  neglected  and  are  the  cause  for  the  secular 
motion  of  tho  spin  .axis.  The  spin  axis  in  Explorer  IV  was  observed  to 
move  as  much  as  three  degrees  per  d.iy  and  to  have  reached  a  maximum 
deviation  from  the  initial  orientation  of  15  degrees  on  the  celestial  spJicre. 
In  .addition  to  the  gr.avity  torque  there  can  .also  be  a  magnetic  torque  acting 
on  the  satellite's  span  .axis.  It  is  saf*!  to  say,  serious  st'cular  movements 
of  tho  spin  axis  of  a  staliilized  s.atellito  do  occur  whicli  will  influence  the 
tirag  aiul  const-quontly  the  di?c,iy  of  the  satellite. 

As  an  instructive  example  jn'vertheless,  consider  a  space  st.ihil- 
iz('d  satellite.  The  importance  of  the  .satellite  orientation  will  then  he 
<lemonHtrat<!<l.  Consider  a  sattjllite  suclt  as  1958  /.eta  which  is  an  Atlas 
missile,  Assiime  this  satellite  to  he  sjjin-stabili zed  with  tlie  spin  vector 
.aligned  with  the  flight  palli  at  tho  initi.al  pc-rigee  point.  Approximately 
one  week  later,  the  ]>erigoe  will  have  regre.ssed  some  58  degrees.  Refer¬ 
ring  to  Sketch  1  below,  it  e.m  be  scon  th.at  tlie  drag  coefficient  will  have  a 
value  somewhere  lietween  the  two  extremes,  A  w.ay  of  representing 
tlie  actual  drag  would  be  to  .assuinc  that  tho  drag  coefficient  depends  upon 
tho  argument  of  the  perigee. 

\  |cos(wo-w)| 
min  \  max  min/ 

For  this  particular  example 

Cj^=1.70+0.65|cos(u)q— w)|  {il.Z4) 
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Sk«iU;h  1 

Tliu  million  of  tho  arj^uiiiont  of  tlio  pcri^co  obviaunly  can  attribute 
a  maximum  error  of  19  per  <;eni  on  the  150th  revolution  whcnug  —  w  a 
90  cU!(<roeM.  'I'liis  hIiouIiI  show  that  the  drag  eoijffieicnt  is  dependent 
upon  tlu!  orbital  elementH  ami  ahould  not  be  included  as  a  constant, 

'I’ci  a  mucli  lesser  degree!,  the  same  can  be  said  for  an  oarth- 
stabiliiiied  satellite.  In  this  case  however,  the  body  axis  will  move  with 
respect  to  the  velocity  victor  only  a  few  degrees  during  the  portion  of 
high  drag  encounter.  'I'liis  example  closely  approximates  the  concept 
of  a  constant  area  and  ilrag  coefficient, 

(c)  Estimation  of  lift  coefficient 

In  addition  to  the  drag,  lift- in  many  instances  will  also  be  present. 
Depending  upon  the  vehicle  orientation  relative  to  the  wind  velocity,  the 
lift  will  act  to  either  increase  or  decrease  the  centripetal  force.  Generally, 
it  will  also  h.ave  a  component  in  the  direction  of  the  side  and  drag  forces. 

The  drag  due  to  lift  will  usually  be  a  small  fraction  of  the  total  drag.  Even 
though  tlic  side  force  due  to  lift  is  also  small,  it  can  create  serious 
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perturbation  in  tlie  orbital  motion.  This  section  will  deal  with  the  mag 
nitude  of  the  lift  coefficient  and  determine  the  effect  it  has  on  the  side 
and  drag  forces. 

The  first  consideration  is  of  a  stabilized  satellite.  The 
orientation  of  the  vehicle  in  this  instance  will  be  fixed  in  apace.  Con¬ 
sider  the  axis  system  of  Sketch  2: 


Sketch  2 


Tlie  satellit<-  orientation  is  in  reference  to  the  r,  f,  or  co- 
ordin.'ites.  For  a  .spin  sl;ibiliz«‘d  sat^rlHte,  the  .angular  momentum  vector, 
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n.,  it)  fixed  in  Bpaco. 

The  momentum  vector  ia  located  with  rcBpect  to  the  N-E-y  inertial 
axis.  Using  Euler  axis  transformation,  the  momentum  vector  can  be 
specified  in  terms  of  the  or-w-n  axis.  The  transformation  is  as  fol¬ 
lows  on  the  next  page.  . 
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Th«!  njomonluin  vector  has  now  been  located  with  I’cspect 
to  the  orbit  pcrif»ee.  Further  transformations  arc  required  to  locate 
it  in  terms  of  the  actual  8atel]it(!  position.  The  satellite  is  locat(!d  in  the 
orbital  piano  by  the  angular  measure  from  the  perigee.  In  this  regard 


li  i  cos  f  II  +  sin  f  II„ 
r  w  n 

•  (<i.25) 

11,  =  cos  f  Il„  —  Bin  f  11 

I  oijo  W 


liic  angular  momentum  vector  has  now  been  located 
with  respect  to  the  satellite  position.  The  velocity  vector  is  also  located 
with  respect  to  the  instantaneous  satellite  position  as 


Vj.  =  V  ^  —  r  e  cos  O'  (2.26) 

V  =  r  e  sin  o  cos  4) 
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J  |h'  iiiiKle  of  attack,  yaw,  and  roll  can  now  be  determined 
as  the  angle  between  momentum  and  velocity  vectors.  The  angle  of  at¬ 
tack  is 


'^f  • 


U,27)a 


Similarly,  the  angles  of  roll  and  yaw  arc 


yaw  =  11 


V  .  H 
r  r 


{^.<17)b 


roll  »  C 


V  .  H 
(Y  a 

V  I  .  Ill  I 


(z.n)c 


lafl  fortes  tNin  be  created  on  unsymnietrlc;  shapeH.  One 
usually  thinks  of  the  flat  plate  or  cone  at  an  angle  of  attack.  Following 
the  same  presiinf ation  as  for  drag,  the  lift  coefficient  for  diffust:  refl ac¬ 
tion  from  a  flat  plat<?  is 


tHJW  A  ,  ,  "Jlir  sin  A  cos  A 

C.j  = — ^2—  erf  {S  sin  A)  + - ^ - 

'fl  '  ‘  r 


The  lift  coefficient  for  the  case  of  specular  reflection  is 


Ij  =  cos  A  j  — ——  sin  A  exp(  —  sin*  A)  + 

's  L  s 


{2.29) 


+  4(sin*  A  f  2^^  crf(S  sinA) 


The  actual  lift  coefficient  lies  between  these  two,  i.c,, 


=  A'Cj^  +  {K  -  DC, 


(Z.IO) 


The  liftine  force  is  always  smaller  than  the  drag  force 


TH  -IbS 


at  moderate  molecular  speed  ratios.  As  an  example,  for  S  s  f»,  the 
L/D  §  0.20,  Though  the  lift  force  is  small  relative  to  the  drag,  its 
effect  in  the  radial  direction  may  be  important  during  the  final  stages 
of  decay, 

2,  Perturbing  functions 

The  preceding  sections  have  served  to  formulate  the 
basic  quantities  needed  to  state  the  force  equations.  It  is  the  purpose 
of  this  section  to  combine  the  various  factors  of  gravity,  drag  and, 
where  applicable,  lift  into  a  form  amenable  for  later  use  in  determining 
the  motion  of  a  near-earth  satellite.  The  equations  will  be  written  with¬ 
out  proof  in  the  orbital  plane  coordinate  system. 


in 


~  +  J  3sin^  <jj  sin*  n)  — 

~T  ^  <1*  O'  +  3)  j 


li 

m 


-jpv 


2JKMK* 


(V^  -  r  cos  n) 


(2,31) 


n 


^ -  Hill  cos  4>  sin^  <r  y  KM  sin^  at  (})  cos  (|i 


—  12  sin* 


a  sin  ({)  cos  (|i 


r  n  e  cos  a )  + 


h 
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F  2JRL  n 

£  = - km  sin  <j)  cos  a  Bin  a  -  Y  [28  sin*  ([)  sin*  a  cos  a 


m 


(C^A\  (  C^N 

12sin  4>  si«  «  cos  o]  —  -  p  si**  a  cos  ""my  I  ^  ) 


III.  EQUATIONS  OF  MOTION  FOR  A  NEAR-EARTH  SATELLITE 


The  force  terms  and  the  perturbing  functions  have  been  developed 
in  the  preceding  section  and  are  now  to  be  equated  to  the  rate  of  change 
of  momentum.  In  gtnieral,  the  cqu.ation8  of  motion  in  spherical  coordi¬ 
nates  are  given  as 


wiiero 


<b  * 

!;  ([)+  0  COB  O' 

o  -  ^  8i»  Ol 
^  r  sin  <t> 


(3.1) 


(3.2) 

(3.3) 


These  equations  are  non-linear  second  order  differential  equations 
and  conceivably  could  be  integrated  numerically  as  they  stand.  This, 
however,  would  be  in  complete  disregard  of  the  brillant  work  of  the  clas¬ 
sical  astronomers.  The  classical  astronomers  developed  several  pertur¬ 
bation  techniques  which  when  applied  to  the  equations  of  motion  offer  many 
particular  advantages.  The  prominent  special  perturbation  techniques 
which  are  applicable  no  the  motion  of  a  near-earth  satellite  are  Cowell's, 
Encke's,  and  variation-of -parameter  methodsiCowell 's  method  is  not  a 
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perturbation  technique  but  is  grouped  as  such  when  comparing  the  relative 
merits  of  these  three  prominent  methods). 

The  choice  between  the  three  methods  depends  upon  the  particular 
trajectory  to  bo  encountered,  Cowell's  method,  which  is  a  direct  integra¬ 
tion  of  the  accelerations,  is  suited  for  large  perturbative  forces.  For 
such  forces  such  as  rocket  thrust,  Cowell's  method  will  afford  the  least 
number  of  computations  per  step  but  will  require  many  steps  (Uef'.  11). 

It  has  the  disadvantage  of  requiring  a  large  number  of  signific.ant  figures. 
Encko's  method  is  best  suited  for  modortate  perturbations  acting  only  during 
;i  segment  of  the  trajectory  such  as  thrus^t  forces  applied  intermittently 
wliile  in  orbit  to  provide  st.abilir.ation  or  in  orbit  transfer  and  turning. 
l'incki-‘'s  tnetliod  is  a  true  perturb.'ition  technique  wherein  only  the  iiccelera- 
tion  difference  between  the  .'ictual  .and  rciference  orbit  is  integriited.  This 
intt.'gratioii  of  only  Mur  dep.arlurt!  from  a  reference  orbit  permits  using 
fewer  signific.ant  figures  .and  l.arger  integration  steps,  Tlu*  natural  im- 
provemmit  to  Encke's  method  was  tlie  use  of  ,a  i:onfinuouHly  elianging 
r<•f('r^.•n(■e  orliif  .  Tliis  is  known  .as  the  v.ariation- of -pa  ramtiler  method  and 
is  suitiMl  for  small  perturh.alion  forcei;  a<!ting  throughout  the  orbit.  This 
tnetliod  permits  large  integrtition  steps  hut  tlu*  eomput.alions  .ari*  I'xten- 
sive  ill  e;u;h  ittep.  As  tliis  study  is  eoiireriied  with  tlu*  satellite's  iiuition 
iu*a  r  .and  ;it  its  demise,  tlie  v.i  riation -of -|>;i  rame  te  r  melltod  affords  tlu* 
minimum  of  computational  lim«  .ami  It.as  flexiliiiily  .and  control  of  llu!  associ- 
.ited  I'omjtutatioiial  errors. 

Tile  va riation -v>f -p.ar.ameter  metliod  is  the  computation  leclinitiue 
used  in  this  study.  'I'lu*  reference  motion  of  tlu*  .Hatellitc  is  represented 
by  a  set  of  orbit  parameters  wliich,  in  tlu*  aljsence  of  the  cirag  forces, 
would  lje  constant  from  rcvululion  to  revolution.  'I'lie  drag  forces,  how¬ 
ever,  cause  these  parameters  to  vary,  and  tlie  differential  equations  of 
motion  ;ire  derived  from  these  parameters.  These  differential  equations 
are  ilicn  integrated  to  describe  the  satellite  motion.  The  description  of 
the  reference  orbit  is  provided  in  the  follov/ing  paragraphs. 

The  reference  orbit,  sometimes  referred  to  as  an  oscillating  or 
instantaiu'ous  ellipse,  is  specified  by  tlirec  parameters.  In  place  of  the 
usual  eccentricity  and  semi -major  axis,  the  dimensionless  parameters 
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p  and  q  are  used.  The  parameter  p  is  the  semi-latus  rectum  of  the  oscil¬ 
lating  ellipse  and  also  represents  the  product  of  velocity  and  radius. 


where 


(roV  )* 


(3.4) 


(3.5) 


Also, 


(3.6) 


The  second  parameter,  q,  is  a  product  of  the  eccentricity  and  semi- 
latus  rectum. 


q  =  P  e 


(3.7) 


One  other  param«!tcr  is  required  to  completely  specify  the  refer¬ 
ence  ellipse.  'I'his  is  the  angle  w  which  locates  the  argument  of  the 
perigee  with  respect  to  some  reference  direction.  The  relation  of  w  to 
the  polar  aiiglc  and  argument  of  the  perigee  is  indicated  l>clow  and  in 
Figure  6. 

w  «  (h  -  f  (3.8) 


where  f  is  the  true  .anomaly,  the  angular  displacement  in  the  orbit  plane 
measured  from  the  perigee. 

The  orientation  of  ti>e  reference  orbit  in  space  remains  to  be 
specified.  This  is  conventionally  done  by  specifying  the  angle  of  inclina¬ 
tion  O',  the  line  of  ascending  node  SI,  and  the  epoch  time.  These  six  or¬ 
bital  elements  completely  specify  the  dynamical  state  of  the  satellite  and 
inherently  require  the  momentum  and  the  energy  associated  in  the  refer¬ 
ence  ellipse  to  be  identical  with  those  of  the  actual  trajectory. 

The  second  order  differential  equations  of  motion  (3.1)  have  an 
equivalent  form  as  six  first  order  differential  equations  in  the  oscillating 
elements.  In  terms  of  these  elements,  the  equations  of  motion  in  the 
oscillating  elements  are  given  (Ref.  13)  as 
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(a) 

dt 

1 

1 

-  r  *  V - Tp 

Ia. 

1  sin  <ti  cos  a 

r  \  nr 

^  % 

(b) 

dP 

...ip 

(c) 

dn 

w 

.sinJ^a) 
AV  '  na  / 

(d) 

da 

a? 

COB  4>  j  ^or  j 

n  ^•"1 

(e) 

dq 

d(i> 

"  COS  f  f 

Si  sin  f 

(f) 

do) 

d4) 

dn 

coa«- 

1  dP  ,  ,  S,  , 

—  -TT-  sin  f - *-  cos  f 

q  d4.  q  _J 

whore  the  quantity  Sj  ia  defined  as 


Those  six  differential  equations  have  been  programmed  for  an 
IHM  70‘i  with  tlie  following  auxiliary  equations. 

KM  r  1  ,  T  / 1  I  j  i  A  j  2  r  1  +  cos  f  +  e*  1 

V^,  =  L  1  W  -pr  (1  -  38in*  ^  sin*  a)J  J 

■  ""V  .  f 

r  =  — s-i —  sin  f 
•M) 

V|^  =  h*  +  (V^  —  rfl  e  cos  a)*  +  (rW  e  sin  a  cos  i}))* 

F'g  _  F  a  1 
m  ~  m  sin  a 
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These  equations  when  programmed  along  with  the  various  sub¬ 
routines  for  evaluating  the  drag  coefficient  and  density  provide  the 
numerical  solution  to  the  equations  of  motion.  The  variation-of-parameter 
method  employed  does  have  the  disadvantage  of  singularities  for  the  case  of 
circular  orbits  (  e  =.0).  The  equations  so  chosen  have  as  the  independent 
variable  the  polar  angle  (j)  which  is  in  effect  the  true  anomaly.  At  the 
terminal  phase  of  the  .satellite  trajectory  where  the  motion  is  nearly  verti¬ 
cal.  this  independent  variable  la  a  poor  choice  since  small  changes  in  it 
are  accompanied  by  large  changes  in  the  dependent  variable.  Ibis  how¬ 
ever  will  occur  below  100  km  and  for  all  purposes  the  satellite  has  ceased 
orbiting.  In  view  of  this,  the  equations  are  not  altered. 


IV.  NUMEUICAl.  SOLUTION 

'i'lic  equations  of  motion  (3.V)  and  the  various  au.Kiliary  equations 
liave  b(!<‘n  j)i’<)graiumt!d  for  an  HIM  709  using  Portram  language.  Ihe 
differential  equatiotis  ;irc  Integr.atod  u>sing  the  Adams -Moulton,  Kungt)-Kutta 
teelinique  (Kef.  1.1).  This  subroutine  integrates  n  set  of  N  simultane¬ 
ous,  first  order  differential  equations.  The  clioice  of  tliis  subroutine 
lies  in  tlu!  available  option  of  int(*gration  stirp- sii'.e.  Tlie  A<lamH- 
Meulton  er ror -coni nil  permits  adopting  a  variable  .step-si/, e,  mode  0, 
or,  a  fixed  step  size,  mode  1,  both  with  the  fourth  order  Kunge-Kulta  in¬ 
tegration  method.  Through  the  use  of  a  variable  step-sine,  self-determined 
by  Ibe  magnitude  of  rouiul-off  error  ami  truncation  error  permitted, 
macliine  time  ean  bo  saved  at  no  loss  of  ciesinid  accuracy.  The  user  is 
free  to  specify  tlic  truncation  error  which  is  equivalent  to  specifying  the 
number  of  significant  figures.  Tlie  round-off  error  control  is  achieved 
willi  the  use  of  double  precision  internally.  The  deriv.Hive  calcuhalions 
are  performed  in  single  j>recision. 

Tlie  variable  step-si.ie  is  limited  in  range  between  a  lower  and 
upper  bound.  The  lower  bound  is  chosen  as  a  compromise  between  the 
desired  accuracy  and  machine  time.  Usually,  the  lower  bound  is  chosen 
100  times  the  upper  bound.  After  considerable  e.xpcrience  is  obtained 
with  the  program,  the  fixed  step-size  mode  1  can  be  used  which  is  five 


TR  “tbS 


33 


ALTiTUOe  (KMJ 


2.02  2.03  2.04  £09  2.0«  207  20* 

0R*0  COEFFICIENT 

Figure  7.  Specular  drag  coefficient  for  a  sphere  as  a  function 
of  altitude 
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times  faster  than  the  variable  mode.  Increments  as  large  as  30  degrees 
have  been  used  with  small  loss  of  accuracy,  as  indicated  in  the  compari¬ 
son  below  (Case  1), 


Pin 
!• 

C 
i: 

K.  A. 
f 

I )  rn  ii 
h 

.  Inc loiiientH  as  large  as  ninety  degreos  definitely  cannot  bo  used. 

Tile  accuracy  cheek  Is  made  with  an  inclined  orbit  witli  no  drag. 
The  luacliine  is  reciuircd  to  repeat  all  pertinent  values  for  twenty  revolu¬ 
tions.  This  wa.'i  achitjvod  through  the  sixth  decimal  place.  Similarly, 
the  rate  of  change  of  right  aseemsion  and  argument  of  perigee  for  an  orbit 
with  no  drag  were  checked  with  available  values  determined  from  series 
solutions  (Kef.  14)  and  agree<l  with  the  desired  aci:uracy  for  tlie  second 
order  oartli  [lotential.  Several  ca.ses  of  results  are  prescmled  below, 
(Case  it). 


EQUATORIAL  ORBIT, 

,0 

m 

f 

h 

1! 

''4- 

u> 

Q 

M 

0 

77  1.8550 

0.09999999 

8,1569862 

0 

0 

0 

1800.000 

72  1.8.519 

0.01000001 

8.1569862 

0.07861)1 

0 

-  0.0392223 

7199.999 

72  1.8557 

0,09999999 

8.1569860 

0.3144134 

0 

-  0.1569119 

POLAR  OKBIT,  —  =  0 
m 

1800.000 

72  1.8703 

0.09999774 

8.1569718 

-  0.0197609 

1.5707963 

-  0.0000000 

POLAR  ORBIT,  —  =  7.5X  lO"'*  km*/Kg 
m 

Ifinp.OOO  ZJV7993  0.C997155  8.15597]!  -  0,0197676  1.570963  -  0.0000000 


A,],  B  r 
rail 

6793.11676  Ion 
0,(11  ''O',.;,! 
l(l3.(t<V'IUiO  Ion 
‘j.HIUKMl  l•.•|l! 

B.!  I.Hl  ilnj. 

0.67.7,11  B  X  10  *  Niiwlonn/Kii. 
1 16.01119  Ion 


A,).  .1  to- 
17. '1318  rad 
6793.877.3  Ian 
0.0  I  B0B77 
10.3.836909  lun 
S.HHH0.3  I  rail 
57  1 .81 

0.6,77I!OX  lO"'  Ni'wloiia/l<|l 
116, OHM  km 
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t  The  program  has  the  following  features,  The  atmospheric  density 

is  obtained  from  a  table  rather  than  an  exponential  equation.  By  specify¬ 
ing  shape,  the  decay  can  be  computed  for  either  cylindrical  or  spherical 
satellites  or  for  the  case  of  constant  drag  coefficient.  For  the  spherical 
or  cylindrical  satellite,  J^he  drag  coefficient  is  computed  on  the  basis  of.. 

.  .  . .  . -  the  molecular  scale  temperature  and  specified  constant  accommodation 

coefficient.  The  drag  coefficient  changes  only  a  fcw_per-ccnt_with-altl-- 
-  tude  as  can  be  observed  in  Figure  7.  These  few  per  cent,  however,  can 
cause  an  additional  height  loss  of  one  kilometer  for  every  ten  revolutions 
as  shown  below  for  the  Case  3  with  A/m  =  IX  lO"^  km*  /Kg. 


o 

II 

2 

N  =  1 0 

Cj^  =  1.8 

C,3  =  2.0 

f 

0 

3604.74 

3604.68 

3604,63 

r  (km) 

6578.1450 

6566.8535 

6  564,8657 

6563.7740 

h  (km) 

aoo.oooo 

188.7886 

186.7408 

185.6494 

c 

.0400000 

0.0140897, 

0.0134895 

0.0148746 

c  (km) 

134.447854 

93.847344 

88.41899 

85,514985 

K.  A. 

0 

-  0.000004 

-  0.000004 

-  0.000005 

Phi (deg) 

0 

64.8345 

64.8 545 

64.8345 

A  much  larger  effect  is  shown  in  the  linear  cccontricily .  For  dif¬ 
fuse  drag  coefficient  the  variation  can  be  greater  deptnuUng  upon  the  defini¬ 
tion  of  the  re-cmission  speed  ratio. 

Normally,  interest  is  shown  only  for  the  values  of  the  orbital  ele¬ 
ments  at  perigee  and  apogee.  However,  at  low  altitudes,  the  perigee  is 
ill  defined  and  it  becomes  necessary  to  define  the  orbital  elementals  at 
every  integration.  The  altitude  for  the  change  in  information  output  can  be 
selected  but  it  has  been  found  to  occur  no  higher  than  l-'O  km  though  it 
may  be  lower. 
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V.  AN  APPROXIMATE  METHOD  OF  SATELLITE 
DECAY  PREDICTIONS 


With  the  aid  of  the  present  compute  program,  ji'^isiona.havc  —  _ 
-been  made  in  the  approximate  satellite  decay  predictions  (Ref .4),  This 
program  has  been  valuable  in  evaluating  certain  functions  and  determin¬ 
ing  relationships  between  orbital  elements.  A  brief  review  will  be  made 
of  the  approximate  method,  including  the  aids  furnished  from.the_computing. . 
program.  A  tost  of  the  approximate  method  will  also  be  presented  on  a 
number  of  earth  satellites.  The  results  are  very  favorable. 

The  basis  of  the  approximate  method  lies  in  the  l.'iw  of  conserva¬ 
tion  of  energy.  The  drag  force  is  assumed  to  occur  instantaneously 
during  a  given  revolution.  The  energy  lost  due  to  the  drag  is  then  sub¬ 
tracted  from  the  satellite  energy,  and  the  orbital  element  adjusted  to 
coiuiicnsate  for  this  loss.  The  loss  of  energy  per  revolution  can  bo  ex¬ 
pressed  as 


dEo 

"dfr 


1  ,  "^1  c 

,3 


21: 


-  r/G 

exp  (If 


(5.1) 


The  approximations  made  in  the  equation  are  within  an  order  of 
e* .  The  equation  employs  an  exponential  density  variation.  The  density- 
scale  height,  G,  is  computed  from  the  1959  ARDG  Model  Atmosphere 
(Uef.6)  and  is  plotted  in  Figure  8.  The  irregularities  between  170  km 
and  280  km  are  due  to  the  variation  in  molecular  temperature. 

Expressing  the  orbit  radius,  r,  as 

r  sqf  c  (1—  co8f)  +  C^(e*) 
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the  integral  can  be  evaluated  in  terms  of  Bessel's  function.  The  loss 
of  energy  per  revolution  is  then  shown  to  be 


where  the  argument  k  represents  the  product  of  density  scale  height,  G, 
and  linear  eccentricity,  c,  the  terms  _ 


»\/~2iTk~  Jo  (k) 

exp(K) 


and 


V  2itTc  1 1  (k) 

cxpliw 


deviate  from  unity  only  as  k  -•0.  For  values  of  k  >4,  the  terms  are 
<1.01.  FQrk>>4,  Equation  (5.2)  is  written  as 


(CnA)KM(l-.)p|/f:  [l 


1  +  0  + 


0.125  -  0.375 


e 


(5,3) 


To  obtain  the  variation  in  the  orliit  duo  to  thi.s  energy  Ions, 
differentiate  the  total  energy  associated  with  the  initial  orbit. 


4^4=+  KMm 
dN  2(q  +  c)* 


(5.4) 


Thus  a  change  in  orbit  energy  will  affect  the  perigee  height  and 
linear  eccentricity.  From  the  two  equations,  the  motion  of  perigee  and 
linear  eccentricity  can  be  evaluated.  Frequently,  the  perigee  and  linear 
eccentricity  are  combined  Into  one  parameter,  the  semi-major  axis 
from  which  the  orbital  period  is  determined. 


q  s  a  —  c 


ZttKM 

—JJT 


(5.5) 


TR  465 


39 


<•'1  ,  'I'  fn I  V^TT^  Jilk 

<1N  '  HN  -  <M  “  V  V'Fn  7  '  '  — F'MT(in“  ~Wn''  FxRi 


od  by 


This  can  be  related  directly  to  rate  of  decrease  in  orbital  peri- 


dP  _  3  P  da  ..  . 

a^r  "  I  a  dN 

The  success  of  orbit  decay  predictions  based  on  the  rate 
of  decrease  in  orliital  period  suggests  that  tlie  two  variables  probably 
arc  independent  for  at  least  large  values  of  linear  eccentricity.  There 
is  no  simple  equation  which  can  be  used  to  separate  dq/DN  and  dc/dN. 
Tlie  r;itio  dq/dc  can  be  approximated  with  varying  orders  of  accuracy. 
The  important  factor  .about  this  ratio  is  its  independeticc  from  the 
satellite  si/c  and  sliapc,  though  it  does  depend  upon  density.  With  this 
in  mind  the  variation  of  dq/dc  with  the  linear  eccentricity  was  computed 
numerically  using  tlie  program  outlined  previously.  The  results  are 
Iilottcd  in  Figure  9  as  1/1  I  (dq/dc)  which  has  limits  of  0  and  1.  This 
information  can  tie  used  directly  to  separiUe  lliu  varialjles.  The  next 
step  in  lorming  an  approximate  equation  for  iiredicling  decay  is 
integration  of  Liquation  (5.6).  Integration  with  limits  extending  from 
the  initial  linear  eccentricity  to  /.ero  or  the  start  of  spiral  decay, 
depending  tipon  initial  altitude  results  in  Equation  (5.8). 

’  m  (.;»[  -  (b.b, 

The  product  (  1  +  (dq/dc))  VcT  is  proportional  to  »/e"  for  all 
values  of  c  >  150.  So  in  fact,  dq  and  dc  are  independent  for  all  values  of 
c  >  160.  As  will  be  shown  later,  values  of  c  <  1  50  are  obtained  only  near 
the  very  end  of  the  satellite  lifetime. 

Tlie  orbital  lifetime  is  obtained  after  furtlier  integration  and  is 
expressed  as 


^  >  1 50 


The  density  parameter,  p  n/C,  for  various  altitudes  is  presented 
in  Figure  10. 
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Figure  11.  Lifetime  of  an  earth  satellite  in  a  circular  orbit 
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Figure  12.  Straight  line  approximation  to  linear  eccentricity  vs  N 


An  independent  exact  solution  is  possible  for  the  asymptotic  case 
of  c  B  0.  For  the  circular  orbit, 


Cj,A 

m 


Na 


(5.10) 


The  solution  is  plotted  in  Figure  1 1  for  a  family  of  load  factors. 

In  this  figure  the  load  factor  of  Sputnik  V  rocket  was  determined,  Tbe 
effect  of  the  earth*  s  oblatcness  on  the  altidude  of  the  satellite  has  an 
important  effect  upon  the  ballistic  coefficient  value  determined. 

A  measure  of  the  validity  of  any  assumption  is  in  the  result.  For 
the  case  of  Equation  (5.8),  it  was  found  in  comparing  the  actual  satellite's 
motion  and  the  machine  computations,  that  the  prediction  became  poorer 
as  the  linear  eccentricity  increased.  It  was  first  thought  that  since  only 
terms  of  order  c  were  used,  the  equation  should  not  be  expected  to  be 
useful  beyond  600  or  800  revolutions.  It  was  found,  however,  than  an 
empirical  factor  of  (c/100)®***‘  could  be  included  in  Equation  (5.8) 
which  then  resulted  in  good  agreement  with  decayed  satellites.  The  em¬ 
pirically  corrected  equation  is  then 


^T  (tm)  '  0  150  (5.11) 

This  equation  has  been  plotted  in  Figure  12.  With  either  Figure  12  or 
Equation  (5.11),  the  satellite  lifetime  and  load  factor  can  be  determined 
using  two  values  of  liacrtr  eccentricity  separated  by  a  known  number  of 
revolutions. 

In  the  absence  of  reliable  satellite  data  at  low  eccentricities, 
machine  computations  were  used  to  evaluate  Equation  (5.11).  Equation 
(5.1  l)was  found  to  predict  too  few  a  number  of  revolutions.  This  could  be 
corrected  by  using  the  (1  +  dq/dc)  factor  which  must  be  corrected  to 
equal  zero  at  c  >  200.  The  final  equation  is 
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This  equation  agreed  with  numerical  calculations  and  satellite  data  at 
all  values  of  linear  eccentricity.  In  Figure  12,  a  comparison  is  made 
of  Equations  (5.11)  and  (5.12).  It  is  seen  they  are  identical  for  c  >  150 
as  expected.  The  nonlinear  term  in  Equation  (5.12)  is  the  source  of 
longer  lifetime  at  low  eccentricities  than  predicted  by  lineariaied  theory. 

Equation  (5,12)  is  also  plotted  in  Figure  14  using  a  larger  scale 
in  order  to  show  the  motion  of  1961  Zeta  and  1959  Epsilon  11.  Those  two 
objects  wore  plotted  .•■ccurding  to  their  observed  altitude  and  linear  ec¬ 
centricity.  In  this  comparison  and  others  which  follow,  tlie  source  of 
information  is  the  Space  Track  Dulletijis.  It  should  bo  recogni/.ecl  that 
this  orl)il  information  is  preliminary  and  in  most  cases  is  revised  after 
more  uhservationa  are  received  and  further  analysis  made.  The  re- 
vised'orhit  information,  however,  usually  lags  far  behind  the  satellite's 
lifetime.  The  working  analysis  sliould,  thtirefure,  he  tlesigned  for  tlte 
preliminary  dat.a. 

'I'Im'  appro-ximate  eciu.ilion  ('>,1/.)  as  fitted  to  satellite  data  is 
ploKi'd  in  h’lgures  13  and  14  as  a  tnnelimi  of  linear  eccentricity.  ‘J'he 
ri'sull  is  a  family  of  curves  f<.»r  each  t)f  the  v.arious  iiltltudcs.  There 
is  a  slight  (lecrease  in  slope  ;it  low  eccentricities  with  increasing  .al¬ 
titudes.  The  noolineai'  term  .added  .at  low  ei'ceiitricilies  is  a  function 
of  attitude.  'J'hc  family  of  curves  can  he  reduced  to  a  single  curve, 
exciiuling  the  nonliiii.'.'ir  term,  u.sing  t/ie  iiorm.ali/.ing  height  f.aclor  ^ 
as  plotted  in  Figure  15.  .Slight  corrections  are  present  for  low  eccentrici- 
tic's. 

Since  it  ha.s  been  jio.ssible  to  collapse  the  family  of  curves  into  one 
curve,  tl'.e  total  number  of  rovobjtions  and  ballistic  coefficient  can  be 
sep.arated.  Using  curves  of  total  number  of  revolutions  versus  linear  ec¬ 
centricity  the  ballistic  coefficient  can  be  determined.  This  has  been  done 
in  Figure  10  for  a  large  number  of  satellites.  The  ballistic  coefficients 
for  each  satellite  is  listed  in  Table  5-1. 
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gurc-  13.  Satellite  lifetime  prediction  as  a  function  of 
altitude  and  linear  eccentricit/ 


Figure  1-}.  Comparison  of  obfiei'ved  and  predicted  satellite 
lifetime  for  various  altitudes 
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Comparison  of  observed  and  predicted  satellite 
lifetime  for  various  Cj^A/w 
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The  satellites  plotted  in  Figure  16  have  inclined  orbits  ranging 
from  28  to  82  degrees.  Their  perigee  heights  range  from  178  to  300 
km.  Even  with  all  these  wide  differences,  the  variation  of  revolutions 
and  linear  eccentricity  were  similar  and  could  be  plotted  on  the  same 
curve.  Figure  16  demonstrates  the  close  approximation  of  Equation 
(S.12|  with  actual  satellites.  The  decay  of  1960  Kappa  appeared  unique 
due  to  its  marked  departure  from  the  empirical  curves.  The  marked 
shift  in  the  data  from  I960  Sigma  should  bo  attributed  to  poor  orbit  in¬ 
formation  at  the  time  of  demise.  The  IBM  program  to  simulate  I960 
Sigma  was  evaluated  at  the  wrong  value  of  area -to -mass  ratio. 


VI.  DISCUSSION  OF  SATELLITE  DECAY 

The  purpose  of  this  study  has  boon  to  develop  techniques  for  the 
analysis  of  satellite  motion  near  its  demise.  It  is  well  known  that  ns 
the  satellite  approaches  its  demise,  the  usual  secular  motions  uiiflorgo 
rapid  elianges  as  a  result  of  the  predomiaiancc  of  drag  force.  The  ap¬ 
proximate  method  developed  and  discussed  previously,  is  sufficient 
in  most  eases  to  predict  the  number  of  revolutions  before  decay.  This 
is  only  a  gross  estimate  and  necessarily  excludes  effects  such  as  the 
oblalt!  earth  and  a  rotating  atmosphere. 

To  explore  in  detail  the  three-dimensional  motion  of  the  satellite, 
tln!  program  just  described  has  been  written.  It  has  not  been  possible  to 
use  this  program  thoroughly  in  an  analysis  and  hopefully  further  use  and 
analysis  can  be  made  at  a  later  date.  In  the  meantime  a  brief  analysis 
was  completed  wliich  permits  a  discussion  of  the  more  prominent 
features  of  the  satellite  motion. 

The  first  orbit  element  to  be  studied  was  the  linear  eccentricity. 
This  would  also  provide  a  check  to  the  approximate  solution.  The  rate 
of  change  of  linear  eccentricity  per  revolution  with  various  ballistic  coef¬ 
ficients  is  presented  in  Table  6-1.  The  variation  of  linear  eccentricity 
for  both  the  sphere  and  the  cylinder  were  also  deterrruned  assuming  no 
heat  transfer  (S^  ±  S)  and  plotted  in  Figure  17.  In  Figure  17,  it  can  be 
seen  that  the  linear  eccentricity  decreases  as  expected  until  the  satellite 
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(lescenda  to  an  altitude  of  220  km,  at  which  time,  the  linear  eccentricity 
experiences  short  but  almost  instantaneous  decreases. 

Further  study  of  satellite  motion  at  very  small  linear  eccentricity 
Indicated  poor  agreement  with  the  approximate  method.  Consequently, 
corrections  were  added  to  the  approximate  curves.  The  numerical  cal¬ 
culations  were  the  only  means  of  obtaining  representative  satellite  motion 
at  small  linear  eccentricity. 

Additional  orbit  elements  wore  studied  and  are  also  reported  in 
Table  6-1.  The  motion  argument  of  perigee  in  the  absence  of  drag  can 
bo  described  by  Equation  (3.9).  At  inclination.^  greater  than  63,34  deg, 
the  perigee  moves  opposite  the  satellite  while  at  lower  inclinations,  the 
motion  is  in  the  direction  of  the  satellite.  For  small  drag  forces 
(D/w  <  2.4978  X  lO"*),  no  effect  duo  to  drag  is  noticed.  However,  for 
large. cir.ig  forces,  radical  (jffects  are  observed.  A  sample  case  is  pre¬ 
sented  in  Figure  18  where  at  .about  1.90  km,  (D/w  «  90),  large  sudden 
variations  in  the  perigee  motion  appear.  The  motion,  however,  reni.ainB 
periodic  until  an  altitude  of  110  km  where  the  perigee  then  begins  to  travel 
;it  the  r.atc  of  phi,  the  polar  angle.  It  is  at  this  same  time  that  the  emis- 
sivity  has  reached  a  maximum  and  begun  to  increase  toward  one.  Further 
detail  can  be  seen  in  Figures  19  and  20.  Those  two  curves  arc  different 
and  tlio  dllferenco  in.ay  ho  duo  to  the  initial  heights.  Similar  differcnc cjs 
have  been  observed  in  satellite  decay.  For  instance,  the  unusual  be¬ 
havior  of  I960  Kappa  as  plotted  in  Figure  16a.  The  usual  behavior  is 
spiral  tieoay  at  a  constant  perigee  clist.mce.  A  possible  explanation  for 
I960  Kappa's  bebavio.  is  assuming  the  ballistic  coefficient  increased 
in  some  manner. 

Tliu  eccentricity  undergoes  unusual  changes  at  low  altitudes.  It 
is  seen  that  the  elliptical  orbit  gradually  shrinks  to  some  magic  mini¬ 
mum,  which  appears  also  to  be  a  fuiictlon  of  initial  conditions,  and  then 
it  rapidly  increases  toward  one  as  the  satellite  path  becomes  ballistic. 

An  example  of  this  behavior  is  presented  in  Figure  21.  Table  6-2  pre¬ 
sents  some  observed  minimum  values  of  eccentricity  during  several 
IDM  runs. 
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Figure  18.  Motion  of  argument  of  perigee  as  a  function  of  polar  angle 


ALTITUDE  (KM) 


TRUE  ANAMOUV  (DEG) 

Figure  19.  Motion  of  true  anamoly  at  low  altitudes 
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gure  23.  Error  in  linear  eccentricity  encountered  by  C_  =  2.0  as  compared  to  cylinder 


Table  6-2 

Minimum  Values  of  Eccentricity 


«i 

^1 

c, 

X 

Cj^A/  w 

E  : 
min 

c 

min 

154 

0.05 

343 

1.3 

0.001309 

8.51 

194 

0.0097 

64.2 

0.87 

0.003825 

21 

169.8 

0,00101 

6.64 

0.87 

0.0007846 

5.1 

200 

0.01483 

99.00 

1.0 

0.00384 

25 

150 

0.00761 

50.054 

1.0 

0.003  • 

19 

200 

0.00075 

4 .966 

1.0 

0.0004455 

2.92615 

2  50 

0.0037676 

2  5 

1.0 

0,00107 

7.028 

2  60 

0.0074871 

50 

1.0 

0.0012 

8.5 

140 

0.044 

300 

1.0 

0.0015 

9.68 

200 

0.00754 

50 

1.0 

0.0027 

18.03 

150 

0.0151 

100 

1.0 

0.001 

9.38 

The  motion  of  the  perigee  and  perigee  altitude  can  be  seen  in 
FiguruH  il  and  23.  The  very  unusual  behavior  of  the  perigee  height  is 
<Ku!  to  tlie  earth's  oblatcncss.  Figure  24  demonstrates  the  rapidity  of 
altiliule  loss  with  large  drag  forces.  This  curve  is  an  over  exaggeration 
of  real  satellites  as  it  is  based  on  free  molecular  flow  even  at  low  al¬ 
titudes  ^)f  100  km.  In  reality,  the  drag  coefficient  is  smaller  than  pre¬ 
dicted  by  free  molecular  theory, 

Tlio  effect  of  drag  on  the  argument  of  right  ascension  is  also 
not  noticeable.  With  larger  drag  forces  there  is  a  very  small  increase 
in  the  backward  motion  of  the  right  ascension.  Even  at  very  large 
forces,  the  change  in  the  motion  of  arguimnt  of  right  ascension  is  small 
and  nearly  unnoticeable. 

In  all  of  the  above  computations,  the  atmosphere  is  assumed  to 
rotate  with  tl\e  earth.  This  wind  velocity  when  added  to  the  satellite 
velocity  produces  an  increased  drag  force.  The  enlarged  drag  force 
augments  the  rate  of  motion  of  all  the  elements.  A  comparison  is  shown 
in  Table  6-1  ,between  the  wind  on  and  wind  off.  The  increase  In  rate  in 
some  cases  is  as  high  as  ten  per  cent.  The  assumed  wind  velocity  has 
not  been  verified  and  quite  possibly  it  is  too  large.  The  assumption 
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becomes  poorer  with  increasing  altitude. 

The  initial  conditions  for  Figures  19  through  Z4  are  identical. 

As  a  final  comparison,  the  approximation  often  made  of  a  con¬ 
stant  drag  coefficient  equal  two  was  compared  to  one  computed.  Figure 
2  5  demonstrates  that  =  2.0  is  a  poor  assumption.  It  resulted  in  a 
loss  of  0.9  km/ rev  linear  eccentricity  and  a  0.3  km/rev  loss  in  altitude 
over  that  for  the  computed  drag  coefficient  using  Equations  (2,18) 
through(2.20}.  These  arc  serious  errors  and  point  the  need  to  use  the 
more  accurate  representations  for  drag  coefficient. 


VII.  CONCLUSIONS 

A  discussion  of  satellite  motion  near  the  end  of  its  lifetime  has 
been  presented.  The  motion  was  observed  to  be  very  irregular  under 
large  drag  forces.  The  effect  of  the  irregular  motion  on  the  total 
satellite  lifetime  was  am, all.  During  the  time  of  the  irregular  motions, 
the  salellite  is  expected  to  have  no  more  than  10  to  30  revolutions  re¬ 
maining  before  decay.  This  is  a  matter  of  one  or  two  days.  For  satel¬ 
lite  decay  predictions  better  than  this,  numerical  calculations  such  as 
those  used  here  should  be  employed.  The  program  written  for  this 
study  is  well  suited  for  the  short  range  (N  <  500  revolutions)  prcdic- 
llotis.  Long  range  predictions  probably  require  too  much  machine  time. 

Considerable  effort  was  made  to  improve  the  approximate  method. 
Satellite  observations  and  numerical  computations  were  used  to  adjust 
the  approximate  method.  The  resulting  equation  when  used  for  dec, ay  pre¬ 
dictions  is  accurate  to  three  significant  figures.  The  equation  graphed 
is  .a  very  good  w.ay  to  handle  the  primary  satellite  orbit  data.  On  the 
graph  ii  is  clear  in  what  manner  the  satellite  decays.  The  significance 
of  the  empirical  factor  (c/100)®’***  required  to  correct  the  approximate 
equation  has  not  been  studied.  It  would  be  interesting  to  attempt  to  dis¬ 
cover  the  effect  it  has  on  the  assumption  made  at  the  beginning  of  the  ap¬ 
proximate  method. 
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